function avg_gdop = AvgGdop(T_range,S)
% in  目标位置范围 T_range [x_l x_r y_l y_h]
% in  监测站位置 S 

xt = T_range(1):1:T_range(2);
[~,size_xt] = size(xt);
yt = T_range(3):1:T_range(4);
[~,size_yt] = size(yt);

[X,Y] = meshgrid(xt,yt);
GDOP = zeros(size_xt,size_yt);

x1 = S(1);y1 = S(2);
x2 = S(3);y2 = S(4);
x3 = S(5);y3 = S(6);
% 观测站高度 10 km
z1 = 10;z2 = 10;z3 = 10;
% 目标高度 0 km
z = 0;

sigma_theta = 1*pi/180;% 方位角测量误差
sigma_fi = 1*pi/180;% 高低角测量误差
sigma_s = 0.01; % 站址误差
% 观测误差的协方差矩阵
dV_dVT = diag([sigma_theta^2 sigma_fi^2 ...
               sigma_theta^2 sigma_fi^2 ...
               sigma_theta^2 sigma_fi^2]);
           
for i = 1:size_yt
    for j = 1:size_xt
        x = X(i,j);
        y = Y(i,j);
        % 站1 
        theta1 = atan((x-x1)/(y-y1));
        r1 = sqrt((x-x1)^2+(y-y1)^2);
        fi1 = atan((z-z1)/r1);
        % 站2
        theta2 = atan((x-x2)/(y-y2));
        r2 = sqrt((x-x2)^2+(y-y2)^2);
        fi2 = atan((z-z2)/r2);
        % 站3
        theta3 = atan((x-x3)/(y-y3));
        r3 = sqrt((x-x3)^2+(y-y3)^2);
        fi3 = atan((z-z3)/r3);
        % 系数矩阵
        F = [(cos(theta1))^2/(y-y1) sin(theta1)*cos(theta1)/(y1-y) 0;
            (x1-x)*cos(fi1)*sin(fi1)/(r1^2) (y1-y)*cos(fi1)*sin(fi1)/(r1^2) (cos(fi1))^2/r1
            (cos(theta2))^2/(y-y2) sin(theta2)*cos(theta2)/(y2-y) 0;
            (x2-x)*cos(fi2)*sin(fi2)/(r2^2) (y2-y)*cos(fi2)*sin(fi2)/(r2^2) (cos(fi2))^2/r2
            (cos(theta3))^2/(y-y3) sin(theta3)*cos(theta3)/(y3-y) 0;
            (x3-x)*cos(fi3)*sin(fi3)/(r3^2) (y3-y)*cos(fi3)*sin(fi3)/(r3^2) (cos(fi3))^2/r3];
        h1 = [(cos(theta1))^2/(y1-y) sin(theta1)*cos(theta1)/(y-y1) 0;
            (x-x1)*cos(fi1)*sin(fi1)/(r1^2) (y-y1)*cos(fi1)*sin(fi1)/(r1^2) -(cos(fi1))^2/r1];
        h2 = [(cos(theta2))^2/(y2-y) sin(theta2)*cos(theta2)/(y-y2) 0;
            (x-x2)*cos(fi2)*sin(fi2)/(r2^2) (y-y2)*cos(fi2)*sin(fi2)/(r2^2) -(cos(fi2))^2/r2];
        h3 = [(cos(theta3))^2/(y3-y) sin(theta3)*cos(theta3)/(y-y3) 0;
            (x-x3)*cos(fi3)*sin(fi3)/(r3^2) (y-y3)*cos(fi3)*sin(fi3)/(r3^2) -(cos(fi3))^2/r3];
        % 布站误差的协方差矩阵
        dXs_dXsT = sigma_s^2.*blkdiag(h1,h2,h3)*blkdiag(h1,h2,h3)';
        C = (F'*F)\F';
        % 定位误差的协方差矩阵
        P = C*(dV_dVT+dXs_dXsT)*C';
        GDOP(i,j) = sqrt(trace(P));
    end
end
sum_gdop = sum(sum(GDOP));
[row,col] = size(GDOP);
avg_gdop = sum_gdop/(row*col);